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The very first stars to form in the universe heralded an end to the cosmic dark 
ages and introduced new physical processes that shaped early cosmic evolu- 
tion. Until now, it was thought that these stars lived short, solitary lives, with 
only one extremely massive star, or possibly a very wide binary system, form- 
ing in each dark matter minihalo. Here we describe numerical simulations 
that show that these stars were, to the contrary, often members of tight multi- 
ple systems. Our results show that the disks that formed around the first young 
stars were unstable to gravitational fragmentation, possibly producing small 
binary and higher-order systems that had separations as small as the distance 
between the Earth and the Sun. 
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The earliest stages of the formation of the first stars in the universe, often termed Population III 
(Pop III), have been well studied (1-3), with current numerical simulations evolving the col- 
lapsing gas from cosmological to protostellar densities (4, 5). Much of the dynamical evolution 
of the gas during this phase is controlled by the formation of molecular hydrogen, the main 
coolant of the gas as it is dragged into the collapsing dark matter minihalos. The amount of H 2 
formed sets the minimum gas temperature in the minihalos at around 200-300K, resulting in 
the first self-gravitating baryonic cores - the initial conditions for primordial star formation - 
having masses of around 1000 times that of the Sun (M ). 

Until recently, it was assumed that each of these cores formed just a single star, because 
no fragmentation was seen in the simulations during the formation of the first protostar. As a 
result, attempts to estimate the final mass of the primordial stars have concentrated on balancing 
the inward accretion of gas from the collapsing core by the radiative feedback from the young 
protostar, with various calculations predicting a final mass in the range 30-300 M (1, 5-7). 

It has been shown (8) that it is possible for the collapsing baryons to break into two distinct 
parts, each evolving independently to form its own star, thereby limiting the mass reservoir 
available for each component. The fragmentation of the gas arises from the chaotic turbulent 
flows that feed the inner regions of the star-forming minihalos. Numerical simulations suggest 
that this may occur in about one-fifth of all cases of primordial star formation (8). However 
this figure is a lower limit, because the simulations were unable to follow the evolution of 
the gas beyond the formation of the initial protostar or binary system. A different study that 
followed the evolution of the gas at later times (9) has shown that it settles into a disk with a 
radius of around 1000 astronomical units (AU), which is unstable to gravitational fragmentation. 
However, the limited mass resolution of this study, and the fact that it did not include the effects 
of the radiative feedback from the new-born stars, rendered its results inconclusive 

Here, we present the results of a high resolution numerical simulation that captures the 
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Figure 1: Density evolution in a 120 AU region around the first protostar, showing the build-up 
of the protostellar disk and its eventual fragmentation. We also see 'wakes' in the low-density 
regions, produced by the previous passage of the spiral arms. 
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formation of the primordial pro to star/disk system from cosmological initial conditions down to 
scales as small as 1.5 AU, and which includes the effects of the accretion luminosity heating as 
the disk builds up around the young protostar. We started by identifying the first dark matter 
minihalo to contain cooling, gravitationally collapsing gas from a simulation of a representative 
cosmological volume (10). 

We then re-zoomed the calculation using a technique called "particle splitting" [employed 
elsewhere in studies of primordial star formation, e.g. (5)] and we focused our attention only 
on the collapsing gas at the center of the minihalo, ignoring the larger-scale evolution of the 
minihalo and its surroundings [see (10) for details]. During this second stage of the simulation, 
the gas collapsed to very high densities. Normally, numerical simulations of this kind stop 
once the gas density exceeds around 10 14 cm' 3 , because the computational cost of evolving the 
entire system beyond this point becomes prohibitively expensive. However, in our simulations, 
we replaced very high density collapsing regions with accreting 'sink' particles (10), each of 
which represents an individual protostar. We then used the measured accretion rate onto the sink 
particles to calculate the luminosity produced by the mass as it fell onto the young protostar. 
This energy was then deposited into the surrounding gas under the assumption that the gas is 
optically thin, thus providing a conservative over-estimate of the heating from the protostar 
[see (10) for details]. 

Prior to the formation of the first protostar, the results of our simulation were very similar 
to those presented elsewhere in the literature (1, 5, 8). However, the use of sink particles allows 
us to follow the evolution of the gas past the point at which the first protostar forms, and hence 
to simulate the build-up of an accretion disk around the protostar (Fig. 1). 

After around 90 years, the disk had nearly doubled in size. For the first 60 years, the struc- 
ture of the disk was dominated by a strong two-arm spiral pattern, a feature common to simu- 
lations of present-day star formation (11). Spiral structures of this kind are a signature of self- 
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Figure 2: Radial profiles of the disk's physical properties, centered on the first protostellar core 
to form. The quantities are mass-weighted and taken from a slice through the midplane of the 
disk. In the lower right-hand plot we show the radial distribution of the disk's Toomre parameter, 
Q = Csk/tiGY,, where c s is the sound speed and k is the epicyclic frequency. Beause our disk 
is Keplerian, we adopted the standard simplification, and replaced k with the orbital frequency. 
The molecular fraction is defined as the number density of hydrogen molecules (% 2 ), divided 
by the number density of hydrogen nuclei (n), such that fully molecular gas has a value of 0.5 
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gravitating disks, in which gravitational torque provides the main source of angular momentum 
transport. Although the spiral pattern started out fairly symmetric, it quickly developed non- 
axisymmetric features. Eventually, the gas in one of the spiral arms became locally unstable and 
started to collapse, forming a second protostar at a distance of roughly 20 AU from the primary. 
The mass of the central protostar at this point was only 0.5 M . 

The surface density remained roughly constant as the disk grew, with the temperature behav- 
ing in a similar fashion and thus the ability of the disk to transport angular momentum remained 
the same as the disk grew. This can be seen by considering the Toomre 'Q' parameter (Fig. [2]), 
which provides a measure of the gravitational instability of the disk. For high Q, the disk is 
stable, while for values around 1, the disk is formally unstable to fragmentation. As the disk 
grew, the value of Q remained around 1 in the outer regions, and so the dynamics of the disk 
were dominated by gravitational instabilities. 

Figure [3} which shows the accretion rate through the disk and envelope as a function of 
radius from the central protostar, helps to explain why the disk became so unstable. The ac- 
cretion rate through the disk was considerably lower than the rate at which material was added 
to the disk. If one assumes a simple a-disk model (12), then the disk accreted with an effec- 
tive a between 0.1 and 1, which is typical for gravitational torques in strongly self-gravitating 
disks (S26). Although this is an efficient mechanism for transporting angular momentum, it was 
nevertheless unable to process the material in the disk quickly enough before more was added 
from the infalling envelope. As a result, the disk grew in mass, became gravitationally unstable, 
and ultimately fragmented. 

For a region of the disk to form a new star, it must be able to rid itself of the heat gener- 
ated as it collapses, and further, it must do so before the motions in the disk shear the region 
apart (S34). In the case of a primordial protostellar disk, this cooling is provided by molecular 
hydrogen. At typical disk number densities, around 10 12 to 10 14 cm -3 , most of the cooling 
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Figure 3: The mass transfer rate through the disk is denoted by the solid black line, while 
the mass infall rate through spherical shells with the specified radius is shown by the dark 
blue dashed line. The latter represents the total amount of material flowing through a given 
radius, and is thus a measure of the material flowing through and onto the disk at each ra- 
dius. Both are shown at the onset of disk fragmentation. In the case of the disk accretion 
we have denoted annuli that are moving towards the protostar with blue dots, and those mov- 
ing away in pink (further details can be found in Section 6 of the online material). The light 
blue dashed lines show the accretion rates expected from an 'alpha' (thin) disk model, where 
M(r) = 3nac s (r) E(r)fT(r), with two global values of alpha and where c s (r), S(r), and 
H(r) are (respectively) the sound speed, surface density and disk thickness at radius r. 
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comes from H 2 line emission. In the centers of the spiral arms, or at number densities of around 
10 14 cm -3 , cooling by collision-induced emission (CIE) from H 2 begins to dominate, and it is 
this process that allowed the initial gravitational fragmentation to take place in the spiral arms 
(see Section 6 of 10). At number densities of around 10 16 cm -3 , the gas becomes optically 
thick to CIE, and begins to heat up further as it is compressed. However, it turns out that pri- 
mordial protostellar accretion disks are unusual in that the dissociation of molecular hydrogen 
plays an important role in regulating their temperature evolution. In contrast to present-day star 
formation, primordial protostellar disks are remarkably hot (Fig. [2]). Any further compression 
causes the hydrogen molecules to dissociate, removing energy from the system. This process 
acts like a thermostat, keeping the temperature in the collapsing region roughly constant while 
the collapse proceeds. Thus, perhaps somewhat counter-intuitively, it is the relatively high tem- 
perature associated with primordial star formation, coupled with a high molecular fraction, that 
allowed the protostellar accretion disk to fragment. 

Fragmentation does not stop with the formation of a second protostar. Only four years after 
the initial fragmentation of the disk, the dense ridge in a neighboring spiral arm also fragmented, 
forming a third protostar. Fifteen years later, the disk fragmented yet again and the protostars 
evolved into a chaotic multiple system. At the stage of the evolution shown in Fig. [T] slightly 
less than 1 M of gas had been converted into protostars, and hence the system was still very 
much in its infancy. Eventually, beyond the epoch followed by our current calculation, the disk 
will lose its ability to fragment once the heating from the protostars becomes strong enough 
to dissociate all of the H 2 , thereby removing the major coolant (15). However the system 
will continue to accrete infalling gas from the collapsing envelope, until one of the protostars 
becomes sufficiently hot to ionize the gas, which finally terminates the accretion (7). 

Are our results representative of primordial star formation in general, or is the gas within our 
halo special in some way? The properties of the gas that assembled the disk in our simulation 
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were remarkably similar to those found in other simulations of primordial star formation (10). 
Crucially, the angular momentum of the inner 4 M - the combined mass of the disk and the 
protostar at the point of fragmentation - was 1.4 x 10~ 3 km s -1 pc, comparable to the 1.2 x 10~ 3 
kms -1 pc reported elsewhere (1). Fragmentation of the kind we have seen should therefore be 
a normal part of the formation of primordial stars, suggesting that the endpoint of the Pop III 
star formation process is significantly more complicated than previously thought. Idealized 
calculations have previously found hints of such widespread binarity among the first stars (16), 
but our work establishes this property for realistic cosmological initial conditions, coupled with 
the radiative feedback, and the extremely high resolution needed to address protostellar disk 
evolution. 

While our simulations do not show how primordial disks evolve beyond the initial disk 
fragmentation, present-day star formation calculations (11,17, 18) predict that as fragmentation 
proceeds in high-mass accretion disks, new protostars form at increasingly larger radii. Because 
gas can only be accreted when its angular momentum matches that of the protostar, it is easier 
for the new objects than for the preexisting ones to gain fresh gas that moves inwards through 
the accretion disk, tending to drive the system towards equal masses (19). If fragmenting Pop III 
systems evolved in a similar fashion, then one plausible outcome of primordial star formation 
would be the formation of nearly equal-mass Pop III binaries. Their potential existence would 
strengthen the case for high-redshift gamma-ray bursts origination from the first stars (20). 
From present-day star formation, we also know that young multiple systems are dynamically 
unstable, over time leading to the dynamical ejection of protostars (21) with the preferentially 
close, high-mass binary systems remaining in the center. If Pop III stars were dynamically 
ejected from such systems before accreting very much gas, then there is the possibility that some 
of these stars may have had masses low enough for them to have survived until the present-day. 
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1 Numerical methods 



The simulations discussed in this paper were performed using a modified version of the Gadget 
2 smoothed particle hydrodynamics (SPH) code (SI). Our modifications involve the addition 
of a treatment of primordial chemistry and cooling, along with a scheme for replacing gas 
in unresolved regions with collisionless "sink" particles (52, S3). We have also implemented 
an approximate model for the effects of the accretion luminosity produced by the accreting 
protostars. Details of these modifications are given below. 

Chemistry and cooling To model the chemical evolution of the metal-free gas, we use a 
chemical network consisting of 45 reactions amongst twelve chemical species: H, H + , H", H^, 
H 2 , He, He + , He ++ , D, D + , HD, and free electrons (S4). For most of these reactions, we use 
the same reaction rate coefficients as in Clark et al. (S4). The exceptions are the three-body 
reaction 

H + H + H^H 2 + H, (1) 
and its inverse reaction, the collisional dissociation of H 2 by H, 

H 2 + H -> H + H + H. (2) 

Considerable uncertainty exists concerning the rate of reaction [T] at low temperatures (56). We 
adopt the rate coefficient given for this reaction by Abel et al. (55), which is the smallest of those 
in common usage, and hence will yield the highest temperature for the dense gas in the center 
of the minihalo (57). This warmer gas is less likely to fragment than the cooler gas yielded 
by alternative choices for the three-body reaction rate coefficient. Moreover, accretion disks 
formed in simulations run with a small value for the three-body rate coefficient are less likely to 
be H 2 -dominated than those formed in simulations using a larger value for the rate coefficient, 
and so again will be more stable against fragmentation. In view of this, we consider our choice 
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of the Abel et al. rate coefficient to be conservative, in the sense that it will tend to reduce the 
likelihood of fragmentation. 

Our decision to use a different rate coefficient for Reaction [T] than in the Clark et al. study 
means that we must also adopt a different rate coefficient for reaction [2} as at high densities the 
ratio between the rate coefficients of reactions 1 and 2 must satisfy 

£ = K, (3) 

where k% and k 2 are the rate coefficients for reactions [T] and EJ respectively, and where K is the 
equilibrium constant for the pair of reactions, given by (S8) 

K = 1.05 X l(T 22 T- ' 515 exp (^ ^ . (4) 

Note that the value of K does not depend on our choice of values for k\ or k 2 , and hence we 
cannot vary these rate coefficients independently: if we change ki, we must also adjust k 2 such 
that Equation [4] remains satisfied. 

Our treatment of the cooling also follows the Clark et al. study. We account for the full 
set of processes described in that paper, including electronic excitation of H, He and He + , 
cooling from the recombination of H + and He + , Compton cooling and bremsstrahlung. How- 
ever, in practice only a few processes play an important role at the densities and temperatures 
characteristic of the protostellar accretion disk - H 2 rotational and vibrational line cooling, H 2 
collision-induced emission (CIE) cooling, H 2 collisional dissociation cooling, and heating due 
to three-body H 2 formation. 

We model H 2 line cooling in the optically thin regime using a comprehensive treatment 
that includes the effects of collisions between H 2 molecules and H and He atoms, other H 2 
molecules, protons and electrons (59), and that accounts for the transition to local thermody- 
namic equilibrium level populations at number densities n ^> 10 4 cm -3 . At densities n ~ 
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10 9 cm -3 and above, the strongest of the H 2 ro-vibrational lines become optically thick, reduc- 
ing the effectiveness of H 2 line cooling. To account for this effect, we use an approach based 
on the Sobolev approximation (S10). We write the H 2 cooling rate as 

A H2 = AKiAii&scui^u, (5) 

u,l 

where n u is the number density of hydrogen molecules in upper energy level u, AE U \ is the 
energy difference between this upper level and a lower level I, A u \ is the spontaneous radiative 
transition rate for transitions between u and /, and /3 CS c,ui is the escape probability associated 
with this transition, i.e. the probability that the emitted photon can escape from the region 
of interest. We fix n u by assuming that the H 2 level populations are in local thermodynamic 
equilibrium, and write the escape probabilities for the various transitions as (S10) 

1 - exp(-r u i) 



Pcsc,ul 

Tul 

where we use the approximation that 



(6) 



r u i ~ a ul L s (7) 

where a u i is the line absorption coefficient and L s is the Sobolev length. In the classical, one- 
dimensional spherically symmetric case, the Sobolev length is given by 

s |df r /dr| ' ^ ^ 

where v th is the thermal velocity, and dt> r /dr is the radial velocity gradient. In our three- 
dimensional simulations, we generalize this as (Sll) 

L s = J^, (9) 
| V • v| 

To prevent the H 2 cooling rate from being reduced by an unphysically large amount in regions 
with small velocity gradients, we limit L s to be less than or equal to the local Jeans length, Lj. 
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We justify this choice of limit by noting that there are strong density gradients in the gas on 
scales L 3> Lj, and so we expect the bulk of the contribution to the H 2 line absorption to come 
from material within only a few Jeans lengths. 

At very high densities (n > 10 14 cm -3 ), CIE cooling from H 2 becomes more effective than 
H 2 line cooling (S12). We model this as discussed in Clark et al., and account for the reduction 
of the CIE cooling rate due to the effects of continuum opacity using the following empirical 
prescription, based on previous high-resolution ID calculations of the formation of primordial 
protostars (S12, S13) 

( 1 — e~ TCIE \ 

A C iE,thick = AciE,thinmin , 1 , (10) 

V r C iE / 



u X 8 . (ID 



where 

J x 10 15 cm" 3 

The reduction in the CIE cooling rate predicted by this prescription agrees to within a factor 
of a few with that measured in a more recent 3D simulation of Pop. Ill star formation (S14), 
but this simulation followed the evolution of the gas only prior to the formation of the initial 
protostar, and hence could not test the validity of the empirical prescription at later times in the 



evolution of the system. In particular, it is unclear whether Eq. 1 1 remains a good approxima- 
tion once material begins to build up in the protostellar accretion disk. To test this, we selected 
a set of SPH particles covering a range of different number densities in the output snapshot pro- 
duced immediately prior to the formation of the second sink particle. We computed continuum 
opacities along forty-eight different lines of sight surrounding each SPH particle by solving the 
equation 

rLi 

Ti = / p(l)Kp(p,T,x H2 )dl (12) 



along each line of sight, where L, is the distance from the particle to the edge of the simulation 
volume, p(l) is the local density at a distance / along the line of sight, T is the temperature at 
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the same point, xu 2 = n^/n, and k p is the Planck mean opacity of the gas. The mean opacity 
is a function of the local density, temperature and chemical composition, and we calculated it 
by interpolation from previously tabulated values (575). An estimate of tcie then follows from 

-y 48 

exp(-r C i E ) = — ^exp(-Ti). (13) 

48 i=l 



We then used Eq. 10 together with this estimate of the CIE opacity to compute the reduction 
in the CIE cooling rate, and compared this with the reduction predicted by our empirical pre- 
scription. The results are plotted in Figure S|4j We see that in practice, the continuum opacity 
begins to reduce the CTE cooling rate somewhat earlier than would be predicted by Eqs. [10] - 
but that at densities above n ~ 10 16 cm" 3 , the two different techniques yield very similar 
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answers. Moreover, the difference between the reduction factor (i.e. the ratio of the optically 



thick and optically thin CIE cooling rates) predicted by Eqs. 10-11 and the one computed from 
our simulation results is never greater than a factor of two over the whole of the regime where 
CIE cooling is important. Given the very strong temperature dependence of the CIE cooling 
(Acie oc T 4 ), a factor of two uncertainty in the cooling rate translates into only a 20% uncer- 
tainty in the gas temperature, which is unlikely to be large enough to significantly affect the 
conclusions of our study. Furthermore, it should also be noted that gas at the typical densities 
of n ~ 10 13 -10 15 cm" 3 found in the accretion disk remains in the optically thin regime as far as 
CIE cooling is concerned. The effects of continuum opacity only become important in gas that 
is already undergoing run-away gravitational collapse. 

Finally, our thermal model also accounts for changes in the internal energy of the gas due to 
the dissociation or formation of H 2 molecules. Every collisional dissociation of an H 2 molecule 
requires 4.48 eV of energy (the binding energy of the H 2 molecule), which is removed from 
the thermal energy of the gas, while every time that a new H 2 molecule is formed by the three- 
body process, the gas gains 4.48 eV of thermal energy. For a gas in chemical equilibrium, the 
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Figure 4: The ratio AciE,thick/AciE,thin, plotted as a function of the number density of H 2 . 



The solid line gives the values predicted by our empirical prescription (Eqs. 10-11), while the 
individual data points show the results of the more sophisticated procedure described in the text. 



simplest way to account for these effects is through the equation of state of the gas. Given 
the hydrogen nuclei number density n and the gas temperature T, the equilibrium molecular 
fraction xn 2 = nu 2 /n can be determined from the Saha equation, and the specific internal 
energy u can then be obtained by summing up the contributions due to the translational motions 
of the hydrogen molecules, hydrogen atoms and helium atoms, the internal degrees of freedom 
of the hydrogen molecules, and the dissociation energy of the H 2 molecules (S16). A convenient 
way to visualize how u varies with temperature is to examine the behaviour of the specific 
heat capacity, u/RT, where R = k^/m-ft, k B is Boltzmann's constant, and m H is the mass 
of a hydrogen atom. In Figure S[5} we show how u/RT varies with temperature for several 
different values of n. At low temperatures, the specific heat capacity is almost independent 
of temperature, but above a temperature of 1500-2000 K, it rises sharply as the H 2 begins to 
dissociate. Physically, this means that once the gas reaches this temperature regime, a much 
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Figure 5: The specific heat capacity u/RT of primordial gas (where u is the specific energy and 
R = k-B/m-fi), computed assuming chemical equilibrium, and plotted as a function of temper- 
ature for three different values of the hydrogen nuclei number density n. At low temperatures, 
the specific heat capacity is almost independent of T, but it increases sharply above T ~ 1500- 
2000 K as the H 2 in the gas begins to dissociate. 



greater amount of energy is required in order to raise the temperature further than would be 
the case in colder gas. We see this effect in our simulations, in the gas undergoing run-away 
gravitational collapse in the disk: despite being strongly heated by compression, it increases its 
temperature only slowly once it reaches the temperature regime in which H 2 dissociation begins 
to occur. 

Unfortunately, although this approach is useful for providing insight into the fhermodynam- 
ical behaviour of the gas, it relies on the gas being in chemical equilibrium, which is not the 
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case for the majority of the gas in our simulations. Chemical equilibrium is reached only for 
gas number densities above n ~ 10 15 cm -3 (S14). Therefore, rather than account for H 2 for- 
mation heating and dissociation cooling by including their effects in our adopted equation of 
state, we instead follow previous studies (S10, SI 4) and include a chemical heating/cooling term 
e c hem = Xh 2 ^h 2 m our equation for the time evolution of the internal energy density of the gas, 
where xh 2 — 4.48 eV and n H2 is the rate of change of the H 2 number density, computed with 
the mass density p held constant. 

Sink particles We create sink particles following the standard SPH prescription (52), and 
adopt a threshold density for sink particle creation nth — 10 17 cnr 3 . At this density, radiative 
cooling of the gas has become almost completely ineffective, and the dominant cooling process 
is collisional dissociation of H 2 (see Section [7] below). We therefore do not expect significant 
fragmentation of the gas to occur at densities n > n t h. Any SPH particle that exceeds our 
threshold density is potentially eligible to be converted into a sink particle, but sink particle 
creation occurs only for particles satisfying several additional criteria. First, the candidate sink 
particle must be located a distance r > 2r acc away from any other sink particle, where r acc is 
the sink particle accretion radius (see below). This criterion prevents sink particle formation 
from occurring within gas that is about to be accreted by another sink particle. Second, we 
check to see whether the smoothing length of the particle is less than the accretion radius of the 
sink particle that it will become. This test ensures that we are resolving the structure of the gas 
within the sink particle accretion radius. Third, we ensure that the candidate sink particle and its 
nearest neighbors are all on the same integration time-step. Once these preliminary criteria are 
met, the dynamical state of the candidate sink particle and its neighbors is examined. In order 
to ensure that sink particle creation occurs only within gas which is undergoing gravitational 
collapse, we require that the total energy of the candidate sink plus the gas within one smoothing 
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length be negative, implying that the collection of SPH particles are gravitationally bound, and 
also verify that the divergence of the particle accelerations is negative. This final check ensures 
that the group of particles is not in the process of being tidally disrupted or bouncing. If all of 
these tests are passed, the candidate particle is converted into a sink particle. Its fifty nearest 
neighbors are removed from the simulation, and their masses and linear momenta are added to 
the sink particle. 

As the simulation progresses, sink particles are allowed to accrete other gas particles that 
move to within r < r acc of the sink particle. However, gas particles within the accretion radius 
are accreted only if they pass several additional tests. First, they must be gravitationally bound 
to the sink particle, and moving towards it. Second, if more than one sink is present, they must 
be bound more strongly to the candidate sink than to any other sink in the simulation. Third, the 
gas particle and the sink must be on the same integration time-step. Once these tests are passed, 
the mass and linear momentum of the SPH particle are added to the sink particle, and the SPH 
particle is removed from the simulation. 

In the simulations presented in this paper, we set r acc = 1.5 AU. Gravitational forces are 
softened on a significantly smaller length scale of 0.066 AU to ensure that the force exerted 
by the sink particle on gas outside of the accretion radius is essentially the same as that of a 
point mass, and to permit sinks to come close enough together to allow us to check for potential 
mergers. For comparison, the minimum SPH smoothing length resolvable in the simulation 
is 0.0551 AU, and so we resolve pressure gradients on a similar scale to gravitational forces, 
as is required to prevent artificial fragmentation (SI 7). We check for mergers between sinks 
by examining whether any two sinks ever come within two protostellar radii of each other, 
taking 50 R Q as a reasonable estimate for the pre-main sequence radii of the Population III stars 
studied in this simulation (see Eq.[T5]below). In this study, we find no evidence for any potential 
mergers. 
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We have also performed a simulation with a larger sink accretion radius of 3 AU. The grav- 
itational softening length in this simulation was again fixed at 0.066 AU. This simulation pro- 
duced qualitatively similar results to those from an r acc = 1.5 AU simulation with the same 
assumed protostellar accretion rate. In the simulation with the larger accretion radius, the first 
signs of fragmentation appeared at a distance of 8 AU from the central protostar, roughly twice 
the separation seen in the simulation with r acc = 1.5 AU. As in the simulations with smaller 
r acc , the disk proceeded to fragment further. The difference between the two simulations is 
easily explained. When the sink particles have smaller accretion radii, we are able to resolve 
fragmentation that occurs on smaller scales. In this particular case, the simulation performed 
using an accretion radius of 1.5 AU was able to model the formation of an object that was swal- 
lowed by the sink particle in the simulation performed with the 3 AU accretion radius. The 
third object to form in the r acc = 1.5 AU simulation is located in roughly the same place as 
the second object in the r acc = 3 AU simulation, demonstrating that the evolution of the disk at 
distances r 3> r acc is largely unaffected by the choice of r acc . 

Accretion luminosity We treat each sink particle formed in our simulation as a separate pro- 
tostar, and account for the energy released by accretion onto the surfaces of these protostars. 
To model the effects of this accretion luminosity, we first compute the bolometric accretion 
luminosity for each accreting protostar 

(7.U..U. 

£acc = * * , (14) 

where M* is the accretion rate onto the protostar, M* is the protostellar mass and is the 
protostellar radius. In this study, we model only the first few hundred years of evolution of 
the protostellar system. Since the Kelvin-Helmholtz relaxation timescale of even a high-mass 
protostar is of the order of a few thousand years (S18), we can be confident that the protostars 
in our simulation will not yet have thermally relaxed, and will still be in the initial adiabatic 
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accretion phase. We therefore relate the protostellar radius to the current protostellar mass and 
mass accretion rate using the following relationship (SI 9) 

(MA 027 / M, \° M 

which was derived for adiabatically accreting, metal-free protostars embedded in a spherically 
symmetric inflow. More recent calculations for the case of rotating infall (S20) find an even 
larger value for during the adiabatic accretion phase (see their Fig. 5), and so using the 
expression given by Equation [15] provides us with a conservative upper limit on the bolometric 
accretion luminosity during this phase. 
Combining Eqs. 14 and [15] we find that 

(MA°- 73 ( \ 059 
L acc ~ 12OOL — - — ^ . (16) 

V M o/ V 10 M ©yr~v 

Using this expression for the bolometric accretion luminosity, we next determine the heating 
rate of the gas surrounding the pro to star from 



where p is the mass density, r is the distance to the protostar, and n-p is the Planck mean opac- 
ity of the gas. We calculate this mean opacity by interpolation, using tabulated values that 
include the effects of both line and continuum absorption, and that account for the extra elec- 
trons provided by ionized lithium (SI 5). Equation [17] assumes that the gas is optically thin to 
the radiation from the accreting protostar. In reality, this is unlikely to be the case. However, 
making this assumption allows us to avoid the extremely high computational cost that would be 
associated with an accurate treatment of the transfer of the protostellar radiation, and also gives 
us a conservative upper limit on the effectiveness of the protostellar feedback. 

To determine L acc , and hence I" 1 *, we need to know the accretion rate onto the protostar, M*. 
Determining this self-consistently during the simulation is difficult, owing to the particle-based 
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nature of our hydrodynamical model. The smallest amount of gas that can be accreted during 
any given hydrodynamical time-step is that corresponding to a single SPH particle, i.e. 1CT 5 M , 
and hence accretion onto the sink particle proceeds as a series of discrete increments, rather than 
a smooth increase in the sink particle mass. To obtain an accurate value for the accretion rate, 
we therefore need to average over a large enough number of time-steps to wash out this artificial 
discreteness. In the simulation discussed in the main paper, we achieve this by recording the 
mass accreted by the sink particle during each hydrodynamical time-step, and determine the 
current value of M* by computing a smoothed average of the accretion rate over the previous 
ten years. For times less than ten years after the formation of the protostar, we average over a 
shorter period, and take the instantaneous accretion rate at the time that the sink particle forms to 
be M* = 5 x 1CT 2 M yr _1 , an overestimate of the true value (see Figure S6b. As the protostellar 
accretion rate typically decreases over time, this procedure gives us a slight overestimate of the 
true rate, as can be seen in Figure S[6} where we compare our estimated accretion rate (denoted 
by the black line) with the actual accretion rate onto the protostar (magenta line), which we 
measure during a post-processing step. In order to explore the sensitivity of our results to our 
method for determining M*, we have also performed two simulations in which the accretion 
rate used within our accretion luminosity calculation is not determined self-consistently during 
the simulation, but is instead maintained at a fixed value throughout. In these two simulations, 
we set M* = 1CT 3 M yr -1 and M* = 10~ 2 M Q yr _1 , respectively, as these values bracket 
the true value of M* during almost all of the period of time simulated, as illustrated in Figure 
S|7} These simulations produce qualitatively and quantitatively similar results to our primary 
simulation, as demonstrated in Section [4] below. 

Finally, we note that in our current simulations, we do not account for the photodissociation 
of H 2 molecules in the surrounding gas by radiation from the central protostar. We justify this 
omission by noting that during the adiabatic accretion phase, the protostar has a very large radius 
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(see Eq.[j~5]above), and hence a low effective temperature. Therefore, only a very small fraction 
of the bolometric accretion luminosity of the protostar is emitted in the form of ultraviolet 
photons energetic enough photodissociate H 2 . As an illustrative example, consider the case of 



an O.5M protostar accreting at arate M* = 1O~ 3 M yr -1 : fromEqs. 15 - 16 above, we find that 
the radius and luminosity of the protostar would be roughly 20 R and 720 L , respectively, 
from which it follows that the effective temperature of the protostar at this point is roughly 
6500 K (where we have approximated the protostellar emission as being from a pure black- 
body). The fraction of the protostellar luminosity that falls within the Lyman- Werner bands of 
H 2 is therefore roughly 10~ 8 , and so the accreting protostar emits roughly 10 39 photons s _1 
within the Lyman- Werner bands. Even if we assume that every single one of these photons is 
absorbed and causes the dissociation of an H 2 molecule, the rate at which H 2 is destroyed by 
photodissociation is less than 10~ 10 M yr _1 , and hence only a negligible amount of H 2 would 
be destroyed during the period covered by our simulations. 

2 Initial conditions 

We begin our study with a cosmological simulation. We model the evolution of dark matter 
in a box with a side length of 200 comoving kpc. The cosmological simulation is initialized 
at z = 99 with a standard ACDM fluctuation power spectrum, with a matter density f2 m = 
1 -Q A = 0.3, baryon density fi b = 0.04, Hubble parameter h = # /100km s" 1 Mpc" 1 = 0.7, 
spectral index n s = 1.0 and normalization as = 0.9 (S21). Following the zoom technique 
used in Greif et al. (522), we first run a coarse dark matter (DM) simulation with 128 3 particles 
and a particle mass of ~ 150 M until the first minihalo exceeding a mass of 5 x 10 5 M 
collapses. We then reinitialize the simulation at z = 99 with additional small-scale power in a 
cube with 50 kpc (comoving) on a side, centered on the location of the minihalo formed in the 
previous step, and replace each DM particle in the original simulation with 512 DM and 512 gas 
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time after formation of first star [yr] 



Figure 6: Evolution with time of the mass of the first sink particle (blue curve), the true accretion 
rate onto the protostar (magenta curve) and the time-averaged estimate of the accretion rate used 
to computing the accretion luminosity (black curve). Our estimate is a slight overestimate of 
the true accretion rate, particularly at early times when the accretion rate is changing rapidly 
with time. 
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Figure 7: Left panel: Evolution with time of the mass of the first sink particle (blue curve) and 
the accretion rate onto that sink particle (magenta curve) in the simulation in which we assume 
an accretion rate of 1CT 3 M yr^ 1 for the purposes of computing the accretion luminosity. The 
true accretion rate remains a factor of a few larger than this assumed rate throughout most of 
the simulation. The rise in the accretion rate visible at t ~ 100 yr is a consequence of the 
same gravitational instability that leads to the formation of the second sink particle. The local 
contraction of the spiral arm causes material to be pulled back in its orbit, towards the collapsing 
region, where it collides with other material. This process causes a sufficient loss of angular 
momentum that the orbit can no longer be maintained and the gas moves inwards, resulting in an 
increased accretion rate onto the central protostar. Right panel: Evolution with time of the sink- 
mass and accretion rate in the simulation in which we assume an accretion rate of 1O~ 2 M yr -1 . 
In this simulation, the true accretion rate onto the sink remains smaller than the assumed value 
throughout most of the simulation. It exceeds the assumed value only at very early times after 
the formation of the first sink, although we caution that the numerical resolution of this initial 
burst of accretion is poor, and that we may be overestimating the rate at which it occurs. In this 
simulation, we also see a rise in the accretion rate associated with the formation of the second 
sink particle, but the magnitude of the increase is much smaller. This is a consequence of the 
fact that the spiral structure close to the protostar is much less pronounced, owing to the greater 
amount of heating produced by the increased accretion luminosity. 
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particles. Outside of this region, a sufficiently large buffer zone is created, where the particle 
mass increases in factors of eight until the mass resolution of the original DM simulation is 
obtained. The final DM and gas particle masses in the high resolution region are ~ 0.3 and 
0.04 M , respectively. We then re -run the cosmological simulation, evolving the gas and dark 
matter until the gas in the minihalo has reached a number density of 10 6 cm~ 3 , by which point 
the gas has gravitationally decoupled from its parent minihalo and has begun to collapse in its 
own right. This occurs at a redshift z ~ 17, at which time the virial mass of the dark matter 
halo is 3.3 x 10 5 M . Note that this value is smaller than the mass quoted above because the 
runaway collapse of the gas cloud occurs slightly earlier than the time at which the minihalo 
was identified in the DM-only simulation. 

At this point we discard the full cosmological simulation and focus our calculation on the 
central collapsing region and its immediate surroundings. These then become the initial con- 
ditions for the second phase of our study. The initial gravitational instability that leads to the 
collapse in the baryons occurs at a number density of around n ~ 10 4 cm~ 3 , in gas that has a 
mean temperature T ~ 270 K, implying that the Jeans mass at this point is around 350 M . 
At the point where we discard the cosmological simulation, the central density is around 10 6 
cm -3 , and so the collapse has already evolved over two orders of magnitude in density from the 
point at which the original instability occurred. To ensure that we capture the entire collapsing 
fragment in our simulations, and to avoid any unphysical boundary effects, we select a spherical 
region containing 1000 M of gas to be re-simulated at higher resolution, and discard the rest of 
the SPH particles. We also discard the dark matter particles from the simulation at this point, as 
the further collapse of the gas is dominated by the self-gravity of the gas, rather than the gravi- 
tational potential of the dark matter, since it contributes only ~ 70 M in the re-refined region. 
As such, removing the dark matter only changes the ratio of rotational to gravitational energy 
by around 7 percent in the region of interest, which is negligible given the variation expected 
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from cosmic variance. 

To account for the effects of the missing gas that should surround this central core, we 
include a external pressure term that modifies the standard gas-pressure contribution to the 
Gadget 2 momentum equation, 



by replacing Pi and Pj with Pi — P cxt and Pj — P cxt respectively, where P cxt is the external 
pressure, and all other quantities have the usual meaning (SI) . The pair- wise nature of the force 
summation over the SPH neighbors ensures that P ext cancels for particles that are surrounded 
by other particles. However at the edge of the cloud, the extra pressure term does not disappear, 
and thus mimics the pressure contribution from a surrounding medium. We use the average 
density and temperature at the edge of our cloud to define the value of P cxt . 

In order to allow us to follow the collapse of the gas up to very high densities, it is necessary 
to substantially increase the resolution. The SPH particle mass in the original cosmological 
simulation was 0.04 M , and so the re-simulated region contains 25, 000 SPH particles. To 
increase the resolution, we 'split' each of these particles into 100 new SPH particles of lower 
mass, giving a new particle mass at the beginning of our re-zoomed simulation of 4 x 10 -4 M . 

The particle splitting is done by randomly placing the sibling particles inside the smoothing 
length of the parent particle (S23, S24, S14). Apart from the mass of the siblings, which is 100 
times less than that of the parent, they inherent the same values for the entropy, velocities and 
chemical abundances as their parents, thereby ensuring conservation of mass, linear momentum 
and energy during the splitting of the particles. Contrary to claims in the literature, splitting 
the particles does not conserve angular momentum by construction, since the mass distribution 
changes during the split (and adaptive mess refinement codes suffer from the same problem). 
However we find the error in the angular momentum introduced during the splitting to be small, 



dvj 
dt 
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around a percent. This process is repeated several times during the simulation to ensure that 
the resolution in the collapsing region remains high. However, in these subsequent 'splits', the 
particle mass is only changed by a factor of five at most. In Table SjTj we list the densities at 
which we split the particles, the temperature and local Jeans mass of the gas at the point where 
splitting occurs, the total mass of gas that undergoes particle splitting at each stage, and the 
particle mass before and after the splitting. We also list the number of SPH particles per Jeans 
mass and per Jeans length prior to the splitting, allowing us to demonstrate that we always apply 
our splitting technique to highly resolved regions. 

Our particle splitting technique ensures that particles of significantly different mass (i.e. 
those with mass ratios greater than 5) never come in contact with one another during the final 
stages of the simulation, minimizing any potential problems due to the mixing of SPH particles 
with wildly different masses. The fact that the inner 82 M of gas is resolved with a uniform 
particle mass of 10~ 5 M also ensures that the disk is formed by SPH particles that are of equal 
mass, and we halt the simulations long before any of the higher-mass particles have reached the 
disk. 

At the point of disk fragmentation, our calculations have in excess of 300,000 SPH particles 
in the disk. Even for mildly self-gravitating disks, this has been shown to be sufficient to 
provide converged evolution for the angular momentum transport (S25, S26). Disks that are 
more strongly self-gravitating, such as the one we present in this Report, are easier to resolve, 
as they have a larger scale height. 

It should also be stressed that rather than promoting artificial fragmentation, the 'kernel' 
softening employed within Gadget-2 to soften gravitational forces on small scales will tend to 
suppress fragmentation in unresolved regions rather than promote it (SI 7, S27, S28), since the 
gravitational forces will be diluted to scales longer than the pressure forces. The high resolution 
that we employ in this study ensures that we capture all of the fragmentation in the density 
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Table 1: We summarize here the details of SPH particle masses that result from the 'particle 
splitting' that is used to increase the spatial resolution in our simulation. The quantities given 
are: the number density at which the resolution is increased; the mean temperature at this 
density; the corresponding Jeans mass; the mass in the region that is refined; the particle mass, 
before and after the split; the number of particles resolving a Jeans mass; and the number 
of particles per Jeans length. For comparison, the bottom two rows show the densities and 
temperatures associated with the onset of gravitational instability in the disk and the formation 
of sink particles. 
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regime in which we have a reliable thermodynamic model for the gas, while in the density 
regime in which our thermodynamical model begins to break down (i.e. n > 10 17 cm -3 ), the 
gas is stable to further fragmentation. 

3 State of the cloud at the onset of star formation 

One potential concern about our findings is that since we have examined only a single realization 
of primordial protostellar collapse, we may be obtaining results that are not typical of Population 
III star formation in general. If the minihalo that we have selected to examine is unusual in some 
way, then the results that we have obtained from it may be misleading. Ultimately, this is a point 
that will need to be addressed by performing a large set of simulations along the lines that we 
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have laid out here, allowing many different realizations to be examined, but at present we do 
not have sufficient computational resources for such a large-scale study. However, as we argue 
in the main article, we can gain confidence that our simulated halo is not atypical by comparing 
its properties to those reported in other studies of Population III star formation. 

As mentioned above, the initial gas distribution for our high resolution SPH simulation was 
taken from a minihalo with a dark matter virial mass of 3.3 x 10 5 M and a dimensionless 
spin parameter A = 0.034. Runaway gravitational collapse of gas within the minihalo occurred 
at a redshift z ~ 17. For comparison, the halos investigated in two of the most recent high- 
resolution studies of Population III star formation (S14, 529) had masses of 5 x 10 5 M and 
5.8 x 10 5 M , respectively, and collapsed at redshifts z = 14 and z = 19, respectively. Our 
halo mass is therefore slightly smaller than these previous values, while our collapse redshift is 
well within the range found in previous work. Our halo spin parameter is slightly smaller than 
the value of A = 0.042 quoted in Turk et al. (529), but lies very close to the mean of the spin 
parameter distribution found in previous studies of high-redshift dark matter halos (S30, S31). 

We can also compare the state of the gas immediately prior to the formation of the first sink 
particle with previously reported results. In Figure S[8j we show plots of the enclosed mass 
M enc and specific angular momentum L of the gas as a function of distance from the centre of 
the collapse (left-hand panels), and of the number density n and radial velocity f ra diai of the gas 
as a function of the enclosed gas mass (right-hand panels). If we compare these radial profiles 
to those reported in previous high-resolution studies (55, S10, S14, 529), we find reasonable 
agreement. We recover the standard p oc r~ 2 2 density profile at radii greater than a few AU, 
and masses above 0.1 M . Our values for the specific angular momentum agree very well with 
the values reported elsewhere (55, 570), demonstrating that our gas is not rotating artificially 
rapidly. We also recover a similar infall velocity profile as in previous studies, with a peak 
value of about 2.5 km s -1 , again in good agreement with previously reported results (55, S14). 
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We do not have sufficient resolution to recover the post-shock, sub-AU scale infall visible in 
Fig. 3B of Yoshida et al. (S14): this would occur within the region that is represented by our 
central sink particle. 

In Figure S(9} we show the radial profile of the three components of the angular momentum 
vector at a point in the cloud's evolution just before the creation of the first sink particle. The 
fact that the relative magnitudes of each component remain constant over such a large distance 
(and enclosed mass) in the cloud shows that the mass shells are, for the most part, rotating 
around a common axis. Only as the enclosed mass approaches a value roughly comparable to 
the original Jeans mass in the minihalo, Mj !init = 350 M , do we see a significant change in the 
orientation of the angular momentum vector. This shows that the collapsing region had a well 
defined rotation axis, while the large-scale motions in the minihalo were much more chaotic. It 
would be interesting to compare these results to those from previous studies, but to the best of 
our knowledge, the required data is not available in the astrophysical literature. 

We have also examined how the product of the mean angular velocity f2 and the local dy- 
namical time tdyn = 1/ \f^Gp evolves as we move away from the site where the first protostar 



forms (Figure S 10). Yoshida et al. (S14) measured this for the central 0.01 M of gas in their 
simulation at the point at which they could no longer follow the evolution of the system, and 
found a value f2td yn = 0.25, in good agreement with our own results. 

Finally, we have also examined the H 2 distribution in the cloud just before the formation of 



the first sink particle. In Figure 1 1 , we show how the H 2 fraction (defined here as the ratio of the 
H 2 number density 7Zh 2 to the number density of hydrogen nuclei, n) varies as a function of n. 
We see that the transition from atomic to molecular hydrogen due to three-body H 2 formation 
occurs at a density of roughly 10 10 cm" 3 , in line with the finding of previous studies that have 
used the same three -body H 2 formation rate coefficient (55, S29). 
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Figure 8: Upper left panel: mass enclosed as a function of the radial distance from the densest 
point in the simulation, plotted just before the formation of the first sink particle. The black 
crosses represent the corresponding values in the simulation presented by Abel et al. (55). Up- 
per right panel: spherically averaged value of the specific angular momentum, L, plotted at the 
same time. Crosses represent the values of L for the indicated enclosed masses in the simulation 
presented by Abel et al. (55). Lower left panel: spherically averaged number density profile at 
the same time, plotted as a function of the enclosed mass. To help guide the eye, the dashed line 
represents the p oc r -2 2 slope found in previous simulations of Population III star formation. 
Lower right panel: spherically averaged radial velocity profile 
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Figure 9: The magnitudes of the individual components of the specific angular momentum 
vector in each shell, plotted as a function of both enclosed mass and radius. The values are 
mass-weighted and are taken from a point in the calculation immediately prior to the formation 
of the first sink particle. The angular momentum is calculated with respect to the densest point 
in the gas. 
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Figure 10: The product of the mean angular velocity f2 and the local dynamical time td yn = 
1 / y/AnGp, plotted as a function of the enclosed gas mass, at a point in the calculation immedi- 
ately prior to the formation of the first sink particle. 
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Figure 1 1 : Ratio of the H 2 number density, nn 2 , to the number density of hydrogen nuclei, n, 
plotted as a function of n, at a time immediately prior to the formation of the first sink particle. 
Note that a value nn 2 /n = 0.5 indicates fully molecular gas. 
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4 Sensitivity to the choice of accretion rate 

Once gas begins to accumulate in an accretion disk around the initial sink particle, its subse- 
quent evolution will depend to some extent upon its temperature, and hence on the degree to 
which it is heated by radiation from the accreting protostar. As already discussed in Section [TJ 
one of the difficulties involved in a self-consistent determination of the accretion luminosity 
of the protostar lies in establishing an appropriately time-averaged estimate for the protostellar 
accretion rate that is insensitive to the particle-based nature of SPH. In our present study, we 
calculate the instantaneous accretion rate onto each protostar by computing a smoothed average 
of the accretion rate over the previous ten years. This procedure yields a reasonable estimate of 
the true accretion rate, as Figure S[6]demonstrates, but nevertheless it remains an approximation. 
It is therefore informative to examine the sensitivity of our results to the value of the protostellar 
accretion rate, in order to help us understand whether the approximation that we are making is 
likely to be a significant source of uncertainty. 

To address this issue, we have examined the results of two additional simulations, performed 
using fixed values of M* = 1CT 2 M yr _1 and M* = 1CT 3 M yr -1 , respectively, for the 
protostellar accretion rate used to compute the accretion luminosity of the protostar. As we 
can see from Figure S|6} these two values bracket the true accretion rate during the period of 
time covered by our simulations. We find evidence for disk fragmentation in both of these 
simulations. When the accretion rate onto the central protostar is low, and hence the accretion 
luminosity is less effective, fragmentation occurs more rapidly than when the accretion rate is 
high. In the M* = 1CT 3 M yr^ 1 simulation, it takes only 105 years before the disk fragments 
to form a second sink particle, whereas in the higher M* simulation it takes 274 years for this to 
happen. Furthermore, when the accretion rate is low, fragmentation occurs closer to the central 



protostar than when the accretion rate is high, as can be seen by comparing Figures S 12 and 
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S 13 (although note the difference in the linear scale of these two figures). Nevertheless, in 
both cases, the general features of the fragmentation are the same: the symmetry in the spiral 
pattern breaks and one arm becomes gravitationally unstable, with part collapsing to form a 
new protostar. The second arm follows quickly after it (no more than a few tens of years in 
both cases), with the result that 4 to 5 stars are formed within a few years of the onset of 
fragmentation in the disk. 

If we compare the results of these simulations with the results of the simulation that used 
a more self-consistent determination of M*, we see that the behaviour of the latter closely 
resembles that found in the M* = 10~ 3 M yr _1 simulation, with fragmentation occurring 
relatively quickly at a distance of roughly 20 AU from the central protostar. This result is 
not particularly surprising, given that the true accretion rate onto the central protostar is only 
slightly larger than 10~ 3 M yr -1 at the point at which the disk fragments (Figure Sj6|, but the 
fact that we still find that fragmentation of the disk occurs even when the adopted accretion rate 
is substantially larger than this argues that any small uncertainties in our determination of M* 
are unlikely to be significantly influencing our results. 

In Figure S fT4| we compare the surface density, temperature, H 2 fraction and Toomre stabil- 
ity parameter for the disks in our three simulations at the point at which the second sink particle 
forms in each simulation. We see that when M* is small, the accretion disk is systematically 
colder at all radii R < 40 AU, and is also significantly denser within the central 5 AU, although 
at larger radii the disk surface densities do not differ by a large factor. The higher central tem- 
perature in the high M* simulation causes a pronounced central dip in the H 2 fraction (defined 
here as n^/n, which means that a value of 0.5 corresponds to fully molecular gas). This dip is 
almost absent in the other two simulations. Comparing the Toomre parameter of the three disks, 
we see that they all have Q ~ 1 at radii greater than about 5-10 AU, indicating that in each case, 
the disk is strongly self-gravitating. It is therefore not surprising that we find fragmentation in 
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Figure 12: The evolution of the density in the M* = 10~ 3 M yr _1 simulation, showing the 
build-up of the accretion disk around the central protostar. Note that the scale and color table 
are different from those used in Figure 1 in the main article. 
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Figure 13: The evolution of the density in the M* = 10~ 2 M yr" 1 simulation, showing the 
build-up of the accretion disk around the central protostar. Note that the scale is different to that 
used in Figure 1 in the main article. 
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Figure 14: Radially averaged disk properties, measured at a point immediately prior to the 
formation of a second sink particle, in simulations with the live accretion rate (black, dot- 
dashed line), M* = 10" 3 M yr" 1 (blue, short-dashed line) and M* = 1(T 2 M yr" 1 (pink, 
long-dashed line). In the low fixed M* simulation, the second sink forms at t$p + 105 years, and 
in the high fixed M* simulation, the second sink forms at tsF + 274 years, where tsF denotes 
the time at which the initial sink particle formed. In the simulation with the live accretion rate, 
the second sink particle forms at t SF + 91 years, a little earlier than in the low fixed accretion 
rate calculation. 
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all three simulations. In fact, the conditions in the disk and infalling envelope are similar to 
those found to lead to fragmentation in a recent study of present-day star formation (S32). 



5 Protostellar masses and accretion rates 



In Figure S 15 we show how the masses of the protostars evolve over the first 120 years of the 
simulation and how the accretion rates onto the protostars vary with time. The first point to note 
is that all of the protostars have a mass of roughly 3 x 1CT 2 M when they form. This is not 
particularly surprising, as this mass scale is simply set by the Jeans mass within the protostellar 
accretion disk. However, accretion from their surroundings rapidly increases the protostellar 
masses, which typically exceed 0.1 M within only 10-20 years of their formation. Another 
point to note is that although accretion onto the central protostar initially proceeds in a fairly 
smooth manner, the onset of fragmentation in the disk produces strong gravitational torques 
that allow more mass to flow into the central protostar, but that also cause the accretion rate to 
become far more variable. Similarly, the complex gravitational interactions between the new 
protostars and the disk also cause the accretion rates onto these protostars to vary strongly with 
time. 

As the disk breaks up, the computational expense involved in following its further evolution 
becomes very large, and so for the present we have been forced to terminate our simulations 
once four or five protostars have formed. At the point at which we stop the simulations, there 
is roughly an order of magnitude difference between the mass of the most massive and least 
massive protostars, but given the high rates at which all of the protostars are accreting gas, plus 
the fact that at this point less than 0.1% of the available gas has been accreted, we cannot state 
with any certainty what the final spread in stellar masses will be, or indeed how many protostars 
will ultimately be formed. 
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Figure 15: Left panel: The evolution with time of the masses of the sink particles. Right panel: 
The accretion rates onto the individual sinks. Once the disk begins to fragment, the strong 
gravitational torques allow more mass to flow onto the central protostar. However, the complex 
interactions between the individual sinks, and between the sinks and the disk cause the accretion 
rates of all of the sinks to become highly variable. 
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6 Rotational support of the accretion disk 

The velocity and mass profiles in the protostellar disk at the point of fragmentation are shown 
in Fig S[T6} We see that despite the strong spiral features, the disk is essentially Keplerian. We 
also see that a substantial fraction of the disk is moving away from the central protostar as the 
new fragment forms, demonstrating that the new protostar is forming in a region that is directly 
responsible for the transport of angular momentum through the disk. In fact, the fragmentation 
occurs at a radius of about 20 AU from the central protostar: roughly the point at which the 
material moving outwards collides with the gas coming in from the envelope. 

The mass profile of the protostar-disk system is also shown in Fig S [T6| We see that the mass 
enclosed by the disk at 20 AU is roughly 2 M , while the mass in the central protostar at this 
time is only around 0.5 M . The disk is therefore significantly more massive than the central 
protostar, a feature commonly reported for the early phases of protostellar mass growth in the 
simulations of present-day star formation (S3 3). 

7 Thermodynamics of the accretion disk 

There are two main criteria that an accretion disk must satisfy before it is able to undergo grav- 
itational fragmentation - the Toomre criterion discussed in the main article, and the Gammie 
criterion, which states that the thermal timescale of the gas in the disk must be a small fraction 
of the orbital timescale (S34) 

i 3 

^-thermal — 3fi = — ^orbital) (19) 

where tthcrmai = e/A, where e is the internal energy density and A is the volumetric cooling 
(or heating) rate, is the rotation frequency of the disk, and £ rbitaJ = 2tt/Q,. If the thermal 
timescale does not satisfy the Gammie criterion, then gravitationally collapsing gas in the disk 
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Figure 16: Upper panel: radial velocity profiles in the plane of the disk, along with the sound 
speed of gas, at a point in the simulation just before the disk fragments to form a second proto- 
star. Regions of the disk that are moving towards the central star are shown by blue dots, while 
those those moving away are shown in red. Lower panel: Mass profile of the protostar-disk 
system just before fragmentation. 
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Figure 17: (a) Dominant heating and cooling processes in the gas that forms the second sink 
particle, (b) Upper line: ratio of the thermal timescale, tthermah to the free-fall timescale, t s , 
for the gas that forms the second sink particle. Periods when the gas is cooling are indicated in 
blue, while periods when the gas is heating are indicated in red. Lower line: ratio of thermal to 
the orbital timescale, t rbitai> for the same set of SPH particles (c) Temperature evolution of the 
gas that forms the second sink (d) Density evolution of the gas that forms the second sink 
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will undergo a thermal bounce and then be sheared apart by the disk rotation, rather than con- 
tinuing to collapse. It is therefore important to understand the thermal evolution of the material 
in the disk, and to verify that the gas which fragments does indeed satisfy the Gammie criterion. 

To do this, we have harnessed the Lagrangian nature of SPH. We first identify the SPH 
particle that finds itself at the center of the region in the disk that collapses to form the second 
protostar. This is the particle that is turned into a sink particle, once it (and its neighbors) pass 
the tests described in SectionfTJ We then go back to a point in the simulation just after the central 
protostar is formed, and follow the evolution of this particle by computing the mean properties 
of the gas as seen by the particle and its 50 nearest neighbors at many different instances in the 
life of the disk. It is this evolution that is shown in Figure S fTTj 

We found that the gas which formed our second sink particle underwent several orbits in the 
disk prior to becoming gravitationally unstable, passing into and out of the spiral density wave 
pattern present in the disk. Gas entering the spiral density wave pattern was compressed, while 
gas exiting it was rarefied, as can be seen from panels 3 and 4 of Figure SjTTJ During periods 
of compression, the dynamical heat input (red circles in panel (a) of Figure S fTTf ) was balanced 
primarily by radiative cooling through H 2 line emission (pale blue line in panel (a) of Figure 
SfTTj), even though the strongest of the H 2 lines all had significant optical depths at this point in 
the calculation. In comparison, cooling from H 2 collision-induced emission (CIE; blue dashed 
line in panel (a) of Figure S fTTj ) made only a minor contribution to the dissipation of energy from 
the disk throughout most of the simulation. The other cooling process included in the Figure, 
H 2 collisional dissociation cooling (dark blue line in panel (a) of Figure SfTTj), is generally 
negligible for the material in the disk. However, the exponential temperature dependence of the 
collisional dissociation rate means that it rapidly becomes more important as the temperature 
rises. 

Once runaway gravitational collapse sets in at t ~ 90 yr, the relative importance of the 
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three main cooling processes changes. As the density and temperature of the gas increase, 
all three processes yield higher volumetric cooling rates. However, the rate of increase in the 
H 2 line cooling rate is slow, owing to the effects of the high line opacity, and so it is quickly 
overtaken by CIE cooling. At the time of the last output dump before the formation of the 
sink particle, H 2 collisional dissociation cooling has less influence than H 2 line cooling or CIE 
cooling, and the gas remains almost entirely molecular, with an atomic hydrogen fraction of 
order 1CT 3 . However, it is plain that the importance of H 2 collisional dissociation cooling is 
increasing rapidly at this point, and we expect it to become the dominant process limiting the 
rate of increase of the temperature of the collapsing gas at higher densities. This is important, as 
this is a non-radiative process, and hence one that will be unaffected by the increasing opacity 
of the collapsing gas. Although we do not follow the details of the collapse beyond our sink 
particle creation density, it is reasonable to assume that the collapse will proceed in a quasi- 
isothermal fashion for as long as there remains H 2 available to dissociate, just as in simulations 
of present-day star formation (S35, S36), or high dynamical range simulations of primordial star 
formation (S14). Only once the H 2 content of the collapsing gas is exhausted will the collapse 
become adiabatic. 



We can also draw an important conclusion from Figure S (T7] regarding the importance of 
accretion luminosity heating in determining the temperature evolution of the gravitationally un- 



stable gas. It is clear from panel (a) of Figure S 17 that the accretion luminosity generally plays 
only a minor role in heating the gas, in comparison to compressional and viscous heating in 
the disk, dominating only when these terms are small or absent. Since our numerical treatment 
of the accretion luminosity heating is designed to maximize its effects, we can be confident of 
finding the same result were we to use a more accurate treatment of the accretion luminosity 
feedback. We can therefore also be confident that radiative feedback from the first star to form 
does not suppress fragmentation in the disk at this stage in the lifetime of the system. We can- 
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not, of course, rule out significant feedback effects later in the lifetime of the disk, when the 
mass of the central sink will be larger, but these lie beyond the scope of our present study. 

As far as the Gammie criterion is concerned, we can see from panel (b) of Figure SfT7|that 
this is satisfied for almost the whole period plotted. Typically, thermal ~ 0.1£ O rbitai, increasing 
above this value only occasionally. Most importantly, thermal Aorbitai becomes small once the 
gas begins to undergo runaway gravitational collapse, decreasing to roughly thermal Aorbitai ~ 
0.01 by the end of the simulation. In contrast, we note that the gas in the disk finds the standard 
Rees & Ostriker criterion for ongoing gravitational collapse (S37), namely that t s > ^thermal* 
more difficult to satisfy. It is this condition that helps to maintain the overall global stability in 
the disk. 

Finally, the results presented here allow us to understand why our conclusions regarding 
the stability of Population III accretion disks differ significantly from those of the previous 



analytical studies {S20, S38, S39). Figure S 17 demonstrates that H 2 line cooling plays a hugely 
important role in the thermal balance of the disk, allowing the disk material to remain relatively 
cold, with a temperature of T ~ 1000-2000 K. However, this process was not included in any 
of these previous analytical studies. They therefore find much higher equilibrium temperatures 
for the gas in the disk. Neglect of H 2 bound-free opacity means that these studies predict inner 
disk temperatures T ~ 6000 K or more, the temperature at which H~ ions first become a major 
source of opacity. At a temperature of 6000 K, the molecular content of the gas is negligible, and 
so the predicted mean molecular weight of the gas in these models also differs by almost a factor 
of two from the value in our cold disks. Together, these effects lead to a significantly higher 
predicted sound-speed for the disk, and hence also a higher Toomre parameter Q. Our simulated 
disks are already marginally stable, and it is likely that a global increase in Q by a factor of a 
few would render them completely stable against fragmentation. The difference between the 
results of these earlier analytical studies and our simulations can therefore be understood as a 
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direct consequence of the difference in the input physics, and underscores the importance of 
properly accounting for the H 2 line cooling. 

8 Potential impact of magnetic fields 

Magnetic fields in the early universe are usually assumed to be extremely weak and their po- 
tential influence on the gas dynamics is neglected in most numerical simulations of Population 
III star formation (S40). Upper limits on the field strength can be derived from the cosmic 
microwave background (S41), from big-bang nucleosynthesis (S42, S43), from modeling the 
reionization history (S44), or from the 21 cm line (S45). Typical values lie around 1CT 9 G in 
comoving units. 

Although a variety of physical processes have been proposed for the creation of a seed mag- 
netic field during inflation, for example via electroweak or QCD phase transitions (S46), most 
studies conclude that the dominant contribution to the magnetic field strength comes from as- 
trophysical processes after recombination. The Biermann battery has been proposed to generate 
fields of the order of 1CT 18 G in the intergalactic medium at a redshift of z 20 (S47). This 
field could be amplified further via the Weibel instability in shocks (S48, S49). The most likely 
mechanism for field amplification, however, is the dynamo process (S50). Indeed, two recent 
studies (S51, S52) conclude that the presence of a small-scale turbulent dynamo is able to in- 
crease the field strength by many orders of magnitude during the collapse of primordial gas 
clouds. Once a disk has formed, the large-scale dynamo and the magnetorotational instability 
may also become important (S3 8, S53). In addition, the presence of coherent fields could launch 
jets and outflows (S54, S55). 

From studies of low-mass star formation in the present day, we know that the presence of 
magnetic fields with field strengths close to the equipartition value can effectively redistribute 
angular momentum via a process called magnetic braking (S56, S57) and can thereby influence 
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the fragmentation behavior (S58, S59). However, this requires the presence of coherent rota- 
tional motions in the disk. Calculations of massive star formation in the solar neighborhood 
with and without radiative feedback (S60, S61, S62) show that the disk around the central proto- 
star quickly fragments to build up a cluster of stars, resulting in extremely complex and chaotic 
gas flows. It is likely that this will reduce the effectiveness with which magnetic fields can 
redistribute angular momentum and drive magnetic tower flows. As a consequence, the frag- 
mentation behavior of the accretion disk around high-mass stars is not changed much by the 
presence of the field. As the disks studied in present-day high-mass star formation calculations 
are similar in nature to the Population III accretion disks modeled in our simulations, we feel 
confident that correctly accounting for the presence of primordial magnetic fields would not 
change our conclusions. 
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